clear
syms A B C D

[input_order]=size(B,2);
[output_order]=size(C,1);

Tp=10;

S=kron(eye(Tp),D);
for i=1:Tp;
    for j=1:(i-1);
        in=(i-1)*output_order+1;
        jn=(j-1)*input_order+1;
        S(in:in+output_order-1,jn:jn+input_order-1)=C*A^(i-j-1)*B;
    end
end

S

all(all(simplify(S.'*S==((S.'*S).'))))